Part 2: Interspecies Competition

# Compute growth ---------------------------------------------------------------
# Load 2008 census data and join additional species data
census_2008 <- bw_census_2008 %>%
  left_join(bw_species, by = "sp")

# Load 2014 census data and filter out resprouts
census_2014 <- bw_census_2014 %>%
  filter(!str_detect(codes, "R"))

bw_growth <-
  # Merge both censuses and compute growth:
  compute_growth(census_2008, census_2014, id = "treeID") %>%
  # Convert data frame to sf object
  st_as_sf(coords = c("gx", "gy")) %>%
  # Clean variables
  select(ID = treeID, species = trait_group, genus, dbh = dbh2, growth, geometry)

SCBI Data

# From forestecology package vignette
# devtools::install_github("rvalavi/blockCV")
library(blockCV)
library(sfheaders)

scbi_2013 <- read.csv("data/scbi.stem2.csv") %>% 
  select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
  mutate(date = lubridate::mdy(date)) %>%
  filter(gx < 300, gy > 300, gy < 600)

scbi_2018 <- read.csv("data/scbi.stem3.csv") %>% 
  select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
  mutate(date = lubridate::mdy(date),
         dbh = as.numeric(dbh)) %>%
  filter(gx < 300, gy > 300, gy < 600)

census_df1 <- scbi_2013 
census_df2 <- scbi_2018
id <- "stemID"

scbi_growth_df <- 
  # Merge both censuses and compute growth:
  compute_growth(census_df1, census_df2, id) %>%
  # Convert growth from cm to mm to make result comperable
  mutate(growth = growth/10,
             sp = as.factor(sp))

# Add spatial information
cv_fold_size <- 100
max_dist <- 7.5

scbi_study_region <- 
 #tibble(x = c(0,400,400,0,0), y = c(0,0,640,640,0)) %>% 
  tibble(x = c(0,300,300,0,0), y = c(300,300,600,600,3)) %>% 
  sfc_polygon()

# Add buffer variable to data frame
scbi_growth_df <- scbi_growth_df %>%
  add_buffer_variable(direction = "in", size = max_dist, region = scbi_study_region)

scbi_cv_grid <- spatialBlock(
  speciesData = scbi_growth_df, theRange = 100, yOffset = 0.9999, k = 9, verbose = FALSE)

# Add foldID to data
scbi_growth_df <- scbi_growth_df %>% 
  mutate(
    foldID = scbi_cv_grid$foldID
  )

# Visualize grid
scbi_cv_grid$plots +
  geom_sf(data = scbi_growth_df, aes(col=factor(foldID)), size = 0.1)

scbi_cv_grid_sf <- scbi_cv_grid$blocks %>%
  st_as_sf()

Create focal vs. competitor data

focal_vs_comp_scbi <- scbi_growth_df %>% 
  create_focal_vs_comp(max_dist, cv_grid_sf = scbi_cv_grid_sf, id = "stemID")

Spread into wide form for modeling (note: should we do this by genus?)

metadata <- read_csv("data/scbi.spptable.csv")

focal_vs_comp_scbi_wide <- focal_vs_comp_scbi %>% 
  group_by(focal_ID, focal_sp, dbh, foldID, growth, comp_sp) %>%
  summarize(comp_basal_area = sum(comp_basal_area)) %>% 
  pivot_wider(names_from = comp_sp, values_from = comp_basal_area, values_fill = 0)

head(focal_vs_comp_scbi_wide)

## # A tibble: 6 x 54
## # Groups:   focal_ID, focal_sp, dbh, foldID, growth [6]
##   focal_ID focal_sp   dbh foldID growth  caca  ceoc  fram  frni  havi  litu
##      <int> <fct>    <dbl>  <int>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1        4 nysy     136.       4  0.103  1.31 0.849 0.204  1.44 0.642  28.3
## 2        5 havi      88        4  0.150  4.58 0.849 0      1.44 2.31    0  
## 3       79 tiam     477.       3 -0.161  4.52 0     3.24   0    0       0  
## 4       80 caca      51.5      4  0.253  2.45 0     0      0    0      13.8
## 5       96 libe      23        3  0.262  1.53 0     4.51   0    0       0  
## 6      101 litu     654.       5  0.552  0    0     0      0    0      31.4
## # … with 43 more variables: nysy <dbl>, saal <dbl>, tiam <dbl>, ulam <dbl>,
## #   cagl <dbl>, unk <dbl>, quru <dbl>, cato <dbl>, fagr <dbl>, libe <dbl>,
## #   astr <dbl>, frpe <dbl>, quve <dbl>, ulru <dbl>, ceca <dbl>, caco <dbl>,
## #   cofl <dbl>, qual <dbl>, juni <dbl>, acru <dbl>, caovl <dbl>, amar <dbl>,
## #   prav <dbl>, vipr <dbl>, prse <dbl>, pist <dbl>, acne <dbl>, qupr <dbl>,
## #   elum <dbl>, juci <dbl>, saca <dbl>, cade <dbl>, crsp <dbl>, ploc <dbl>,
## #   chvi <dbl>, ilve <dbl>, divi <dbl>, crpr <dbl>, qumi <dbl>, qufa <dbl>,
## #   quco <dbl>, rops <dbl>, pivi <dbl>

# write_csv(focal_vs_comp_scbi_wide, "data/focal_vs_comp_scbi_wide.csv")

Plot all trees in SCBI plot:

ggplot(data = scbi_growth_df) +
  geom_sf(aes(col = sp), size = 0.5) +
  labs(title = "All trees in SCBI site")

Plot all the trees (??)

ggplot(data = bw_growth %>% filter(species == "maple")) +
  geom_sf(aes(col = species), size = 0.5)+
  labs(title = "All maple trees in Michigan Big Woods site")

ggplot(data = bw_growth %>% filter(species == "evergreen")) +
  geom_sf(aes(col = species), size = 0.5)+
  labs(title = "All maple trees in Michigan Big Woods site")

ggplot(data = bw_growth %>% filter(species == "oak")) +
  geom_sf(aes(col = species), size = 0.5)+
  labs(title = "All maple trees in Michigan Big Woods site")

ggplot(data = bw_growth %>% filter(species == "short_tree")) +
  geom_sf(aes(col = species), size = 0.5)+
  labs(title = "All maple trees in Michigan Big Woods site")

ggplot(data = bw_growth %>% filter(species == "shrub")) +
  geom_sf(aes(col = species), size = 0.5)+
  labs(title = "All maple trees in Michigan Big Woods site")

maples <- read_csv("https://rudeboybert.github.io/SDS390/static/data/maples.csv")

Model loos like this: \[ \begin{eqnarray*} y &=& \beta_0 + \beta_{\text{dbh}}\cdot \text{dbh} + \beta_{\text{evergreen_bm}}\cdot \text{evergreen_bm} + \beta_{\text{maple_bm}}\cdot \text{maple_bm} + \beta_{\text{misc_bm}}\cdot \text{misc_bm}\\ && + \beta_{\text{oak}}\cdot \text{oak} + \beta_{\text{short_tree}}\cdot \text{short_tree} + \beta_{\text{shrub}}\cdot \text{shrub} + \epsilon \end{eqnarray*} \]

where

\(\text{dbh}\) is the “starting” DBH in 2008
\(\text{bm}\) is “biomass” as estimated by the basal area.

Show the models

model_maple <- lm(growth ~ dbh + evergreen + maple + misc + oak + short_tree + shrub, data = maples)

model_maple can be written as: \[ \begin{eqnarray*} y &=& .120157747 + .008636484\cdot \text{dbh} -0.691935999\cdot \text{evergreen_bm} + .079657822\cdot \text{maple_bm} -2.402421803\cdot \text{misc_bm}\\ && -0.115779874\cdot \text{oak} + .648948573\cdot \text{short_tree} + 2.506065233\cdot \text{shrub} + \epsilon \end{eqnarray*} \] Residual analysis Look at model significance

What factors influence growth? (maybe later)

After you’ve selected the model, describe the relationship between \(y\) = growth and \(\text{dbh}\) and \(\text{biomass}\), showing all your work:

In `model2`:

\[ \begin{eqnarray*} y &=& .120157747 + .008636484\cdot \text{dbh} -0.691935999\cdot \text{evergreen_bm} + .079657822\cdot \text{maple_bm} -2.402421803\cdot \text{misc_bm}\\ && -0.115779874\cdot \text{oak} + .648948573\cdot \text{short_tree} + 2.506065233\cdot \text{shrub} + \epsilon \end{eqnarray*} \]